use "...\itus_diary_F1.dta", clear

foreach var of varlist sleep nap study work_hh out_leisure ind_leisure {

foreach num of numlist 1/24 {
areg `var' sunset_time i.woy if hour == `num', absorb(dist_season) cluster(dist_week)
gen `var'_coef`num' = _b[sunset_time]
gen `var'_cil`num' = _b[sunset_time] - 1.96*_se[sunset_time]
gen `var'_ciu`num' = _b[sunset_time] + 1.96*_se[sunset_time]

}
}

foreach var of varlist sleep nap study work_hh out_leisure ind_leisure {
gen `var'_coef = .
gen `var'_cil = .
gen `var'_ciu = .
}

foreach var of varlist sleep nap study work_hh out_leisure ind_leisure  {
foreach num of numlist 1/24 {
replace `var'_coef = `var'_coef`num' if hour == `num'
replace `var'_cil = `var'_cil`num' if hour == `num'
replace `var'_ciu = `var'_ciu`num' if hour == `num'
}
}


drop sleep_coef1- ind_leisure_ciu24


gen hour_t = .
local x = 1
foreach num of numlist 19 20 21 22 23 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 {
replace hour_t = `x' if hour == `num'
local x = `x' + 1
}

drop hour

ren hour_t hour

*Panel a
twoway (line sleep_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick) ) ///
(rarea sleep_cil sleep_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Sleep: 95% CI") label(1 "Sleep: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )  
		

*Panel b
twoway (line nap_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick) ) ///
(rarea nap_cil nap_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Nap: 95% CI") label(1 "Nap: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )   

*Panel c
twoway (line ind_leisure_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick) ) ///
(rarea ind_leisure_cil ind_leisure_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Indoor Leisure: 95% CI") label(1 "Indoor Leisure: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )   

*Panel d
twoway (line study_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick)  ) ///
(rarea study_cil study_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Study: 95% CI") label(1 "Study: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )   
	
*Panel e
twoway (line out_leisure_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick) ) ///
(rarea out_leisure_cil out_leisure_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Outdoor Leisure: 95% CI") label(1 "Outdoor Leisure: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )     
	
*Panel f
twoway (line work_hh_coef hour,  sort msymbol(O) lpattern(solid) lcolor(black) lwidth(vthick) ) ///
(rarea work_hh_cil work_hh_ciu hour,  sort lwidth(none)  color(black%25)  ) ///
,yline(0, lcolor(black%10)) xline(6, lcolor(black%10)) text(25 14 "Average Daily Sunset Time = approx. 6.20 pm", color(black)) ///
text(-20 3 "Day T", color(black)) text(-20 9 "Day T+1", color(black)) graphregion(color(white)) plotregion(icolor(white)) ///
ytitle("Effect of an Hour Delay in Sunset", axis(1))  xtitle(Hour of Day) legend(order (1 2) ///
label(2 "Work: 95% CI") label(1 "Work: Coef.")) ylabel(30(10)-30,angle (30))  xlabel(1 " " 2 "8pm" 3 " " 4 " " 5 "11pm" 6 " " ///
7 " " 8 "2am" 9 " " 10 " " 11 "5am" 12 " " 13 " " 14 "8am" 15 " " 16 " " 17 "11am" 18 " " ///
19 " " 20 "2pm" 21 " " 22 " " 23 "5pm" 24 " " )     
	